Nonequilibrium magnetic and superconducting phases in the correlated electron 
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A theory is presented for a nonequilibrium phase transition in the two-dimensional Hubbard model 
coupled to electrodes. Nonequilibrium magnetic and superconducting phase diagram is determined 
by the Keldysh method, where the electron correlation is treated in the fluctuation exchange ap- 
proximation. The nonequilibrium distribution function in the presence of electron correlation is 
evoked to capture a general feature in the phase diagram. 
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Introduction — While our understanding of the physics 
of electron correlation has matured, there are still in- 
triguing avenues that are yet to be fully explored. One 
such avenue is, in our view, strongly correlated electron 
systems in nonequilibrium situations. While there are a 
body of intense studies on nonequilibrium states in strong 
AC fields such as strong light sources that can trigger 
photo-induced insulator-to-metal transitions (see [l| and 
refs therein) , or nonequilibrium states in strong DC elec- 
tric fields that can introduce pair-creation of electron and 
holes in dielectric breakdown 0, Q , here we pursue yet 
another situation, where nonequilibrium states are con- 
ceived for an open, correlated electron system coupled to 
electrodes (Figl2] (a) inset). Two effects are expected to 
arise from the bias voltage V across the electrodes. One is 
bi-carrier doping, i.e., electrons and holes are simultane- 
ously doped, since two Fermi energies exist due to the two 
electrodes. Naively one might guess that this can make 
the system superconducting with Cooper pairs formed by 
electrons or holes at half-filling, but this has to be tested. 
There is in fact the second effect of the electron-electron 
scattering in nonequilibrium that makes the originally 
sharp Fermi edges to be smeared. The smearing is ex- 
pected to degrade magnetic orders Q , which in our case 
implies that the smearing should act to reduce antiferro- 
magnetic order. The natural question then is: will this 
also destroy the d-wavc superconducting state? 

So wc study this problem, which is motivated by two 
recent experimental developments. One is the fabrica- 
tion of functional structures with oxides [1, 0, 0|. In 
refs. d, Q, properties such as superconducting transi- 
tion in a clean electron gas formed at an interface of two 
insulating oxides was studied, while Ueno et al. have suc- 
ceeded in controlling the superconducting transition in an 
electrolyte-SrTiOs system by changing the applied volt- 
age. Nonlinear transport properties near the Mott tran- 
sition at interfaces have also been theoretically studied 
in [1, [^, [l3| . The second is an experimental observation 
by Pothier et al. of the nonequilibrium electron distribu- 
tion — double-step Fermi distribution — in a mesoscopic 
copper wire attached to two electrodes [ll| ■ They showed 
that the step in the distribution is rounded due to elec- 
tron scattering. The smearing effect is expected to be 



even stronger in correlated electron systems, so that in 
order to examine the nature of the nonequilibrium phase 
transitions in correlated systems, it is theoretically essen- 
tial to develop a method to deal with the nonequilibrium 
distribution of quasi-particlcs in a self-consistent manner. 
Here we perform this by using the Keldysh method, while 
the interaction is treated within the fluctuation exchange 
approximation (FLEX) [T^ . The transition to a super- 
conducting state is studied with the linearized Eliashberg 
equation. 

We briefly comment on the past studies on supercon- 
ductivity transition out of equilibrium. In a pioneering 
work by Chang and Scalapino [l^l who have solved the 
electron-phonon model self-consistently, it was pointed 
out that nonequilibrium conditions such as irradiation of 
light can cause the quasiparticle distribution function to 
deform, and, under certain conditions, can lead to higher 
Tc as observed in conventional s-wave superconductors 
[3, [ill- niore recent attempts, the critical properties 
near an insulator-superconductor transition was studied 
in [iBl followed by several authors [l3, ■ 

Here we adopt the Hubbard model, a prototype in the 
study of magnetism, superconductivity and other phase 
transitions in correlated electron systems. In the two- 
dimensional square lattice near half-filling, the ground- 
state is the Mott insulator with an antiferromagnetic or- 
der [l^ . When chemically doped with carriers (electrons 
or holes), it is believed that Cooper pairs are formed with 
d-wave symmetry and the system turns superconducting, 
12 .. 20, as also discussed phcnomcnologically in 

23, 0] ■ So the question here is the effect of nonequilib- 
rium on these. 

FLEX+Keldysh method — We consider a thin layer 
of strongly correlated material described by the two- 
dimensional Hubbard model which is coupled to elec- 
trodes. Here we have assumed for simplicity the top 
and bottom electrodes (Fig[2](a) inset), since wc want to 
single out the effect of nonequilibrium situation caused 
by the bias voltages, while a lateral attachment of the 
electrodes would cause a change in the spatial sym- 
metry of the phases. The total Hamiltonian is then 
given by if = TJsys + i?sys-ioad + Hicad, where iJ^ys is 
the Hubbard Hamiltonian with the hopping integral t 
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(taken to be the unit of energy hereafter) and the re- 
pulsive interaction /7, iJsys_iead is the system-electrode 
coupling, and -fficad the electrode Hamiltonian. The 
self-energy is then a sum = Sf^^d + ^int of contri- 
butions from the electrodes and from the interaction, 
where a = r,a,<,>,K denote, respectively, the re- 
tarded, advanced, lesser, greater, and Keldysh compo- 
nents (see e.g., [23). If we label the top and bottom 
electrodes by 7 = 1,2 the electrode self-energy becomes 

YK _ 9,- ^'1 f o,-|"h "-A't yr _ _„■ £t 
'^Icad — ^' Z^7=l,2 2 '"^^^^^ 2T ' ^load — ' Z^7 = l,2 2 

where is the coup ling strength between the system 
and the electrodes, [J, and the energy dependence in 
the density of states is neglected. Here, /i-y represents 
the respective chemical potential of the electrodes while 
their temperature T is kept to be the same, and we fix 
= 0.001. We note that if the coupling is too strong 
(F-y 0.1), no ordering takes place. Within the FLEX, 
the self-energy arising from the electron interaction is 
given by 

s>;<(P,^) (1) 

= / ^^E^cff'<(fc,^')G>'<(p-fc,c.-c.'), 
k 

where p, k are momenta, uj the frequency, N the num- 
ber of fc-points considered. The retarded component 
is determined from ImE*" = ^ (S^ — E^), to which 
the real part is related via Kramers-Kronig's relation. 
This relation between the lesser, greater components and 
the retarded component is valid for other quantities as 
well. The fluctuation interaction is given by = 

C/^Im (|x>'< + ^Xc'^ - Xo '^), where XsiXc) represent 
the spin (charge) susceptibilities, whose retarded compo- 
nents are = x5/(l -Uxh), Xc = x5/(l + Uxh)- Here 
Xo is the irreducible susceptibility, 

Xo-^iQ,^) (2) 
f ^^EG<'>(fc,^')G>'<(fc + 9,c. + c.'). 

The lesser and greater components of y" are deter- 
mined with the help of the Langreth rules [23|. Finally, 
the Green's function is determined from the self-energy 
through the Dyson equation, (C'")"^ = {Gq"')^^ — 
Y^r,a ^Yie retarded and advanced components, and 
G>'< = G"'E>^<G°, for the Keldysh component [H with 
{Gq°')~^ = uj — Ef^ zt iS. The process is repeated until a 
self-consistent solution is obtained. The nonequilibrium 
distribution function fcs can be extracted psj from the 
relation, 

G^ = (l-2/eff)(G^-G"). (3) 

We seek for a self-consistent solution of the above equa- 
tions with iteration until the self-energy converges. In 
the calculation we take a 64 x 64 grid for the square 
Brillouin zone, while an almost logarithmic mesh [23.[29j 
with 301 points for the w-axis is used. We note that the 



distribution function /off deviates significantly from its 
non-interacting form (double step Fermi function), 

/Off = [Fi/fd(c^ - + F2/fd(^ - A^2)]/(Fi + F2), (4) 

(with /fd being the Fermi-Dirac distribution) as an effect 
of strong interaction as discussed below. 

The superconducting transition is studied in terms 
of the linearized Eliashbcrg equation, here extended to 
nonequilibrium. To this end, we iteratively (i = 1, 2, . . .) 
obtain a series of anomalous self-energy {(pf) and anoma- 
lous Green's function (Ff) using E", Xs.c obtained in 
the previous step. With a random initial guess for (pl, 
the Green's function is determined from the linearized 
Nambu-Gor'kov equation, F[ = (jfJ[{LoZf - {ej^+X)\ 
with LuZ = Lo~ [E'~(w) - (E'~(-w))*]/2 and X ^ [S''(w) + 
(E''(— a;))*]/2. Then the lesser (greater) components are 
calculated with the generalized distribution function, 

= {1-2U){F: -Ft). (5) 

We assume here that the distribution of the anomalous 
component is the same as the normal component. Fi- 
nally, we plug this into the Eliashbcrg equation, 

(6) 

k 

where the effective interaction in the singlet chan- 
nel is P>„< ^ uHrn{lx^ < ~\x>'<). The eigen- 
value of the linearized Eliashbcrg equation is ob- 
tained as A = linii^oo ||'/>r+il|/||0nL where \\4>\\\ = 

dujj^^pcjf {p,U})\'^Y^'^ is the norm. The supercon- 
ducting transition takes place when A exceeds unity. 

Nonequilibrium phase transition — We have applied 
the above formalism to obtain the nonequilibrium phase 
diagram of the two-dimensional (square lattice) Hubbard 
model attached to two electrodes by numerically solving 
the equations self-consistently. In equilibrium the phase 
diagram within FLEX as obtained in [l^] has an anti- 
ferromagnetic phase when the doping level 5 = 1 — n is 
small, which is taken over by a d-wavc superconductor 
as 5 is increased. So the interest is how these are modi- 
fied. The result for the nonequilibrium situation in Fig. [1] 
(a), where we plot the spin susceptibility lm.Xs{<iT^) for 
y = 0.1 and a doping level 5 = 0.14, shows that the anti- 
ferromagnetic fluctuation remains strong near half-filling. 
We have four incommensurate peaks around q = (7r,7r), 
as in equilibrium. The effect of increased bias is that 
the peak height is reduced and the peak position on en- 
ergy axis shifts upwards as displayed in Fig. [T](c), where 
ImxP°^'^(q, ijj) for q = (tt, I.Itt) is plotted. We notice that 
no features such as dip or hump appear around oj ^ V . 
The dominant superconducting solution in Fig. [1] (b) 
is again similar to the equilibrium case, that is, the d- 
wave gap has the largest Ad for the linearized Eliashbcrg 
equation. However, the critical temperature Tc at which 
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FIG. 1: (Color online) Spin susceptibility linxsiQ,^^) (a) and 
superconducting gap function Re(/>(fc,t<; = 0) (b) are color- 
coded versus momentum for a bias = 0.1 above the critical 
value, with the doping level S — O.U, U — 4.5, and fi = —0.35. 
Dashed lines in (b) represent the nodes, (c) The peak value 
Imxr''''(q,w) versus a; for V = 0.100 - 0.200 from top to 
bottom for q = (tt, l.lvr). T = 0.002 for (a)-(c). (d) The 
temperature dependence of the Eliashberg eigenvalue Ad for 
the c?-wave pairing for V — 0.00 — 0.10 from top to bottom. 



Xd — I docs depend on V , as shown by the temperature 
dependence of Ad plotted in Fig. [1] (d). So the bias V re- 
duces Tc, until finally the superconducting state no longer 
exists even at zero temperature when the bias becomes 
too strong. We define this as the critical bias Vc- For the 
region of the band filling for which the antiferromagnetic 
order dominates over the superconducting state, we can 
similarly define the bias-dependent Ncel temperature 
as the temperature at which the spin susceptibility di- 
verges d^l • The spin susceptibility is reduced as the bias 
in increased, until antiferromagnetic order ceases to exist 
even at zero temperature beyond the "critical Neel bias" 
Vat. The doping dependence of the Neel bias and the 
critical temperatures for a fixed bias is shown in Fig. [2] 
(a). We can see that, while the antiferromagnetic (AF) 
phase is relatively persistent, the superconducting (SC) 
region rapidly shrinks with the bias V and disappears 
at 1/ ~ 0.1. In Fig. [2] (b) wc give the zero-temperature 
phase diagram on the (V, 6) plane. The Neel bias, peaked 
at the undoped point with Vn — 0.36, decreases with 
the doping, and the AF phase is replaced with the SC 
phase around 5 ~ 0.1 with a maximum critical bias for 
SC Vc — 0.1. As we further increase the doping, the 
SC phase finally disappears. Figure [5] (c) schematically 
summarizes these phase transitions in the (T, V, S) space. 
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FIG. 2: (Color online) (a) The phase diagrams for U/t = 4.5 
for various values of the bias voltage V with AF (antiferro- 
magnetic) and SC (superconducting) phases. Origins of the 
three panels are shifted. Inset: Schematic sample (shaded) 
configuration with two electrodes, (b) The zero-temperature 
phase diagram on the {V, S) plane, (c) Schematic phase dia- 
gram in the (T, V, 5) space. 

Nonequilibrium distribution function — As was exper- 
imentally found in a tunneling measurement in a meso- 
scopic wire of copper by Pothier et aZ.,[ll[ the nonequi- 
librium electron distribution becomes smeared from the 
simple double-step Fermi distribution f^g due to elec- 
tron scattering. In correlated materials, with a strong 
electron-electron interaction, we expect a greater smear- 
ing effect to take place. Indeed, as we shall reveal below, 
the key feature to understand the nonequilibrium phase 
diagram for the open Hubbard model may be captured by 
the way in which the nonequilibrium distribution func- 
tion is rounded by the interaction effect. 

Figure [3] (a) plots the effective distribution /off de- 
fined in eq.(l3]) obtained sclf-consistently for V = 0.06 and 
V = 0.38. The temperature in the electrodes and thus in 
/°ff is set to zero. If we compare the result with the 
corresponding noninteracting distribution function f^g 
(eq.dH)) (dashed lines), fes is seen to significantly de- 
viate from f^g. More importantly, we find here that the 
effective temperature approximation breaks down, that is, 
we cannot fit /off to f^g with the temperature as a fitting 
parameter. Instead, the best fit to the data is given by 

r l-ae(-+^/2)/-, uj<-V/2 
fcs = I -(1 - '2a)Lj/V + 1/2, ~V/2 <uj<V/2 (7) 
[ ae-(--^/2)/-, V/2 < u 

where a and r are the fitting parameters. The parameter 
T having the dimension of energy represents the extent 
to which the distribution is smeared from the double- 
step function. Wc have found in Fig. [3] (b) that the 
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FIG. 3: (Color online) (a) Nonequilibrium distribution func- 
tion for two values of the bias V at half filling (5 = 0). Dashed 
lines are the noninteracting distribution function /^ff. (b) 
Nonequilibrium distribution function (dots) against u in the 
Lo < -V/2 region for V = 0.08, 0.19, 0.32, 0.47, 0.63, 0.80 
from the top, where curves represent a fit with eq.((7ll. (c) The 
smearing parameter r against the bias V for various values 
of 5 and fixed U = 4.5 and T = 0. Fitting errors are smaller 
than the size of each symbol. 



fitting function eq. ([7]) is general in the present open Hub- 
bard model in that all the data with different parameters 
(y, r, U,S, . . .) are reproduced within the numerical er- 
rors. If we specifically plot the bias-dependence of the 
smearing parameter in Fig. [3] (c), we can see that they 
fall upon an approximately linear, universal relation 

T(xV. (8) 

The proportionality constant in this relation depends on 



the interaction strength U and the coupling F to the elec- 
trodes, but not on the filling S as seen from the figure. 
The constant is reduced when the coupling to the elec- 
trode becomes stronger. 

From the viewpoint of the smeared distribution, we 
can view the bias-driven phase transitions in the follow- 
ing way. We have seen in Fig. [2] (b) that the AF (SC) 
orders die out at y ~ 0.4 {V ~ 0.1) respectively. In terms 
of the relation ([5]), these values correspond to the smear- 
ing parameters r ~ 0.1 (t ~ 0.02). We can then note 
that these values are similar to the highest Neel (critical) 
temperatures in the zero bias phase diagram (Fig. [2] (a), 
upper panel). Thus, the transition takes place when the 
smearing parameter r attains a value (depth of each or- 
dered region in Fig. [2] (c) translated to r) similar to the 
transition temperature (height of the region). AF spin 
fluctuations are suppressed in finite bias voltages in this 
manner, which is similar to what happens in itinerant 
electron magnets [1] . 

Discussion — We have obtained a nonequilibrium 
phase diagram for the two-dimensional Hubbard model 
and pointed out the possibility of controlling the phases 
in strongly correlated heterostructures by external bias. 
Both of AF and SC regions shrink with the bias V, which 
we attribute to the smearing of the nonequilibrium dis- 
tribution function. While the smearing can be reduced if 
we make the system more strongly coupled to the elec- 
trodes (in, e.g., a thinner sample), this will lead to the 
destruction of order because a larger electrode coupling F 
will make the spin fluctuations weaker. Thus we conclude 
the smearing of the distribution function is an important 
property of correlated electron systems out of equilib- 
rium, and an experimental verification of this should be 
interesting. We have to make a causion that FLEX em- 
ployed here has limitations in that it ignores the vertex 
correction, cannot address, due to its weak-coupling na- 
ture, the behavior close to the Mott insulator point. Ef- 
fects of electrodes (on, e.g., the pairing symmetry) when 
they are attached laterally are also intriguing. A more 
ambitious future problem is a possibility of bi-carrier in- 
duced superconductivity in nonequilibrium, for which the 
present formalism may serve as a starting point. 

TO wishes to thank Thomas Dahm and Yoichi Yanase 
for helpful advices. 
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